module mMisc

        use mPrecision

        implicit none

        ! Number to parallelize cumulative sum. DEREK---optimize this.
        integer, parameter :: CUMULATIVE_SUM_N_PARALLEL = 16

        interface normalizedGaussian
                module procedure normalizedGaussian_full
                module procedure normalizedGaussian_ll_ul
        end interface normalizedGaussian

        interface arithmeticSeries
                module procedure arithmeticSeries_4
                module procedure arithmeticSeries_8
                module procedure arithmeticSeries_i
        end interface arithmeticSeries

        interface arithmeticPartition
                module procedure arithmeticPartition_4
                module procedure arithmeticPartition_8
        end interface arithmeticPartition

        interface arithmeticVector
                module procedure arithmeticVector_r
                module procedure arithmeticVector_i
        end interface arithmeticVector

        interface cumulativeSum
                module procedure cumulativeSum_i
                module procedure cumulativeSum_r
        end interface cumulativeSum

        interface outerDifference
                module procedure outerDifference_r
                module procedure outerDifference_i
        end interface outerDifference

        interface outerProduct
                module procedure outerProduct_4
                module procedure outerProduct_8
                module procedure outerProduct_i
        end interface outerProduct

        interface upperTriangular
                module procedure upperTriangular_0
                module procedure upperTriangular_up
        end interface upperTriangular

        interface lowerTriangular
                module procedure lowerTriangular_0
                module procedure lowerTriangular_up
        end interface lowerTriangular

        interface iminloc
                module procedure iminloc_4
                module procedure iminloc_8_1D
                module procedure iminloc_8_2D
!                module procedure iminloc_16
        end interface iminloc

        interface imaxloc
                module procedure imaxloc_4
                module procedure imaxloc_8
!                module procedure imaxloc_16
        end interface imaxloc

        interface allocate
                module procedure allocate_r_1D_ptr
        end interface allocate

        interface deallocate_ptr
                module procedure deallocate_r_1D_ptr
        end interface deallocate_ptr

        interface copyToLowerTriangle ! REMOVE
                module procedure copyToLowerTriangleDeprecated
        end interface

        interface copyToUpperTriangle ! REMOVE
                module procedure copyToUpperTriangleDeprecated
        end interface

#ifdef PGF90
        interface erf
                real function erf(x)
                        real, intent(in) :: x
                end function erf

                double precision function derf(x)
                        double precision, intent(in) :: x
                end function derf
        end interface
#endif ! PGF90

        private
        public :: assert
        public :: difference
        public :: normalizedGaussian
        public :: iminloc,imaxloc
        public :: arithmeticSeries, arithmeticPartition, arithmeticVector
        public :: geometricSeries, geometricPartition
        public :: cumulativeSum
        public :: outerDifference, outerProduct
        public :: upperTriangular, lowerTriangular
        public :: copyToUpperTriangleDeprecated, copyToLowerTriangleDeprecated
        public :: copyToLowerTriangle ! REMOVE
        public :: copyToUpperTriangle ! REMOVE
        public :: logBin
        public :: odd, even
        public :: true
        public :: checkAlloc, allocate, deallocate_ptr
        public :: almostEqual

contains

        subroutine assert(logic, errorMessage)

                logical, intent(in) :: logic
                character(len = *) :: errorMessage

                if(.not. logic) then
                        write(0, *) 'ERROR: ', errorMessage
                        stop
                end if

        end subroutine assert

        real(D) function difference(a, b)

                real(D), intent(in) :: a, b

                difference = a - b

        end function difference

        pure real(D) function normalizedGaussian_full(x, mean, rms_variance, norm)

                use mConstants

                real(D), intent(in) :: x, mean, rms_variance, norm

                real(D) :: root2Sigma, toSquare

                root2Sigma = sqrt(2.0_D) * rms_variance
                toSquare = (x - mean) / root2Sigma

                normalizedGaussian_full = norm &
                        / ( sqrt(PI) * root2Sigma ) &
                        * exp( -toSquare * toSquare)

        end function normalizedGaussian_full

        real(D) function normalizedGaussian_ll_ul(x, mean, rms_variance, &
                                                       norm, normLL, normUL)

                use mConstants

                real(D), intent(in) :: x, mean, rms_variance, &
                                       norm, normLL, normUL

                real(D) :: root2Sigma, toSquare

                root2Sigma = sqrt(2.0_D) * rms_variance
                toSquare = (x - mean) / root2Sigma

                normalizedGaussian_ll_ul = norm * 2.0_D & 
                        / ( sqrt(PI) * root2Sigma &
                            * (erf( (normUL - mean) / root2Sigma ) &
                             - erf( (normLL - mean) / root2Sigma )) ) &
                        * exp( -toSquare * toSquare)

        end function normalizedGaussian_ll_ul


        ! Fortran 90, Ch. 22.2.
        pure subroutine arithmeticSeries_4(vector, first, factor)

                real(4), intent(out) :: vector(:)
                real(4), intent(in) :: first, factor

                integer, parameter :: SUBVECTOR_FACTOR = 32 ! DEREK---optimize this one

                ! I found SUBVECTOR_FACTOR = 32 to be optimum for
                ! geometricSeries and
                ! size(vector) = 100 on dual Xeon with ifc -tpp7 -xW
                ! But I also found that ifc -O3 -tpp7 -xW gave a factor of 3
                ! speed improvement with no subvectorizing, i.e., with
                ! THRESHOLD > 100:

                ! When size(vector) > THRESHOLD, subvector scaling will be used.
                ! I'm guessing at the optimum value for this because it's not
                ! so important for me right now:

                integer, parameter :: THRESHOLD = 2 * SUBVECTOR_FACTOR

                real(4) :: subvectorFactor
                integer :: i, il, iu

                vector(1) = first

                if( size(vector) <= THRESHOLD ) then
                        do i = 2, size(vector)
                                vector(i) = vector(i - 1) + factor
                        end do
                else
                        do i = 2, SUBVECTOR_FACTOR
                                vector(i) = vector(i - 1) + factor
                        end do
                        subvectorFactor = factor * SUBVECTOR_FACTOR
                        il = SUBVECTOR_FACTOR
                        do
                                iu = 2 * il
                                if ( iu >= size(vector) ) then
                                        vector( il + 1 : size(vector) ) = &
                                            vector(1 : size(vector) - il) &
                                            + subvectorFactor
                                        exit
                                else
                                        vector( il + 1 : iu ) = &
                                            vector(1 : il) + subvectorFactor
                                end if
                                subvectorFactor = subvectorFactor &
                                                + subvectorFactor
                                il = iu
                        end do
                end if
        end subroutine arithmeticSeries_4

        pure subroutine arithmeticSeries_8(vector, first, factor)

                real(8), intent(out) :: vector(:)
                real(8), intent(in) :: first, factor

                integer, parameter :: SUBVECTOR_FACTOR = 32 ! DEREK---optimize this one

                ! I found SUBVECTOR_FACTOR = 32 to be optimum for
                ! geometricSeries and
                ! size(vector) = 100 on dual Xeon with ifc -tpp7 -xW
                ! But I also found that ifc -O3 -tpp7 -xW gave a factor of 3
                ! speed improvement with no subvectorizing, i.e., with
                ! THRESHOLD > 100:

                ! When size(vector) > THRESHOLD, subvector scaling will be used.
                ! I'm guessing at the optimum value for this because it's not
                ! so important for me right now:

                integer, parameter :: THRESHOLD = 2 * SUBVECTOR_FACTOR

                real(8) :: subvectorFactor
                integer :: i, il, iu

                vector(1) = first

                if( size(vector) <= THRESHOLD ) then
                        do i = 2, size(vector)
                                vector(i) = vector(i - 1) + factor
                        end do
                else
                        do i = 2, SUBVECTOR_FACTOR
                                vector(i) = vector(i - 1) + factor
                        end do
                        subvectorFactor = factor * SUBVECTOR_FACTOR
                        il = SUBVECTOR_FACTOR
                        do
                                iu = 2 * il
                                if ( iu >= size(vector) ) then
                                        vector( il + 1 : size(vector) ) = &
                                            vector(1 : size(vector) - il) &
                                            + subvectorFactor
                                        exit
                                else
                                        vector( il + 1 : iu ) = &
                                            vector(1 : il) + subvectorFactor
                                end if
                                subvectorFactor = subvectorFactor &
                                                + subvectorFactor
                                il = iu
                        end do
                end if
        end subroutine arithmeticSeries_8

        pure subroutine arithmeticSeries_i(vector, first, factor)

                integer, intent(out) :: vector(:)
                integer, intent(in) :: first, factor

                integer, parameter :: SUBVECTOR_FACTOR = 32 ! DEREK---optimize this one

                ! I found SUBVECTOR_FACTOR = 32 to be optimum for
                ! geometricSeries and
                ! size(vector) = 100 on dual Xeon with ifc -tpp7 -xW
                ! But I also found that ifc -O3 -tpp7 -xW gave a factor of 3
                ! speed improvement with no subvectorizing, i.e., with
                ! THRESHOLD > 100:

                ! When size(vector) > THRESHOLD, subvector scaling will be used.
                ! I'm guessing at the optimum value for this because it's not
                ! so important for me right now:

                integer, parameter :: THRESHOLD = 2 * SUBVECTOR_FACTOR

                integer :: subvectorFactor
                integer :: i, il, iu

                vector(1) = first

                if( size(vector) <= THRESHOLD ) then
                        do i = 2, size(vector)
                                vector(i) = vector(i - 1) + factor
                        end do
                else
                        do i = 2, SUBVECTOR_FACTOR
                                vector(i) = vector(i - 1) + factor
                        end do
                        subvectorFactor = factor * SUBVECTOR_FACTOR
                        il = SUBVECTOR_FACTOR
                        do
                                iu = 2 * il
                                if ( iu >= size(vector) ) then
                                        vector( il + 1 : size(vector) ) = &
                                            vector(1 : size(vector) - il) &
                                            + subvectorFactor
                                        exit
                                else
                                        vector( il + 1 : iu ) = &
                                            vector(1 : il) + subvectorFactor
                                end if
                                subvectorFactor = subvectorFactor &
                                                + subvectorFactor
                                il = iu
                        end do
                end if
        end subroutine arithmeticSeries_i

        pure subroutine arithmeticPartition_4(vector, first, last)

                real(4), intent(out) :: vector(:)
                real(4), intent(in) :: first, last

                call arithmeticSeries_4(vector, first, &
                        (last - first) / ( size(vector) - 1 ) )

        end subroutine arithmeticPartition_4

        pure subroutine arithmeticPartition_8(vector, first, last)

                real(8), intent(out) :: vector(:)
                real(8), intent(in) :: first, last

                call arithmeticSeries_8(vector, first, &
                        (last - first) / ( size(vector) - 1 ) )

        end subroutine arithmeticPartition_8

        pure function arithmeticVector_r(first, add, n)
                real(D), intent(in) :: first, add
                integer, intent(in) :: n
                real(D) :: arithmeticVector_r(n)
                call arithmeticSeries(arithmeticVector_r, first, add)
        end function arithmeticVector_r

        pure function arithmeticVector_i(first, add, n)
                integer, intent(in) :: first, add
                integer, intent(in) :: n
                integer :: arithmeticVector_i(n)
                call arithmeticSeries(arithmeticVector_i, first, add)
        end function arithmeticVector_i

        ! Uses subvector scaling technique, as descibed in Numerical Recipes in
        ! Fortran 90, Ch. 22.2.
        pure subroutine geometricSeries(vector, first, factor)

                real(D), intent(out) :: vector(:)
                real(D), intent(in) :: first, factor

                integer, parameter :: SUBVECTOR_FACTOR = 32

                ! I found SUBVECTOR_FACTOR = 32 to be optimum for
                ! size(vector) = 100 on dual Xeon with ifc -tpp7 -xW
                ! But I also found that ifc -O3 -tpp7 -xW gave a factor of 3
                ! speed improvement with no subvectorizing, i.e., with
                ! THRESHOLD > 100:

                ! When size(vector) > THRESHOLD, subvector scaling will be used.
                ! I'm guessing at the optimum value for this because it's not
                ! so important for me right now:

                integer, parameter :: THRESHOLD = 2 * SUBVECTOR_FACTOR

                real(D) :: subvectorFactor
                integer :: i, il, iu

!                factor = exp( log(last / first) / ( size(vector) - 1 ) )
                vector(1) = first

                if( size(vector) <= THRESHOLD ) then
                        do i = 2, size(vector)
                                vector(i) = vector(i - 1) * factor
                        end do
                else
                        do i = 2, SUBVECTOR_FACTOR
                                vector(i) = vector(i - 1) * factor
                        end do
                        subvectorFactor = factor ** SUBVECTOR_FACTOR
                        il = SUBVECTOR_FACTOR
                        do
                                iu = 2 * il
                                if ( iu >= size(vector) ) then
                                        vector( il + 1 : size(vector) ) = &
                                            vector(1 : size(vector) - il) &
                                            * subvectorFactor
                                        exit
                                else
                                        vector( il + 1 : iu ) = &
                                            vector(1 : il) * subvectorFactor
                                end if
                                subvectorFactor = subvectorFactor &
                                                * subvectorFactor
                                il = iu
                        end do
                end if
        end subroutine geometricSeries

        pure subroutine geometricPartition(vector, first, last)

                real(D), intent(out) :: vector(:)
                real(D), intent(in) :: first, last

                call geometricSeries(vector, first, &
                        exp( log(last / first) / ( size(vector) - 1 ) ))

        end subroutine geometricPartition

        ! Computes the cumulative sum of an array, with an optional
        ! scalar added to every element
        pure recursive function cumulativeSum_i(array, scalar) result(ans)
                integer, intent(in) :: array(:)
                integer, intent(in), optional :: scalar
                integer :: ans( size(array) )

                integer :: add
                integer :: i

                if( size(array) == 0 ) return

                add = 0
                if( present(scalar) ) add = scalar

                ans(1) = array(1) + add
                if( size(array) < CUMULATIVE_SUM_N_PARALLEL ) then
                        do i = 2, size(array)
                                ans(i) = ans(i - 1) + array(i)
                        end do
                else
                        ans( 2:size(array):2 ) = &
                                cumulativeSum_i( array(2:size(array):2 ) &
                                + array( 1:size(array) - 1:2 ), add )
                        ans( 3:size(array):2 ) = ans( 2:size(array) - 1:2 ) &
                                + array( 3:size(array):2 )
                end if

        end function cumulativeSum_i

        pure recursive function cumulativeSum_r(array, scalar) result(ans)
                real(D), intent(in) :: array(:)
                real(D), intent(in), optional :: scalar
                real(D) :: ans( size(array) )

                real(D) :: add
                integer :: i

                if( size(array) == 0 ) return

                add = 0.0_D
                if( present(scalar) ) add = scalar

                ans(1) = array(1) + add
                if( size(array) < CUMULATIVE_SUM_N_PARALLEL ) then
                        do i = 2, size(array)
                                ans(i) = ans(i - 1) + array(i)
                        end do
                else
                        ans( 2:size(array):2 ) = &
                                cumulativeSum_r( array(2:size(array):2 ) &
                                + array( 1:size(array) - 1:2 ), add )
                        ans( 3:size(array):2 ) = ans( 2:size(array):2 ) &
                                + array( 3:size(array):2 )
                end if

        end function cumulativeSum_r

        pure subroutine logBin(xs, n, ll, ul, deltaOverX)

                real(D), intent(out) :: xs(:)
                integer, intent(out) :: n
                real(D), intent(in) :: ll, ul
                real(D), intent(in) :: deltaOverX
                real(D) :: onePlusDeltaOverX

                onePlusDeltaOverX = 1.0_D + deltaOverX

                n = ceiling( log (ul / ll) / log (onePlusDeltaOverX) )

                call geometricSeries( xs(1 : n), ll * sqrt(onePlusDeltaOverX), &
                                      onePlusDeltaOverX )

        end subroutine logBin

        ! outerOperations are defined as (A @ B)(i, j) = A(j) @ B(i), for whatever operation @.
        pure function outerDifference_r(a, b)
                real(D), intent(in) :: a(:), b(:)
                real(D) :: outerDifference_r( size(b), size(a) )
!REMOVE                outerDifference_r = spread( a, 2, size(b) ) &
!                                  - spread( b, 1, size(a) )
                outerDifference_r = spread( a, 1, size(b) ) &
                                  - spread( b, 2, size(a) )
        end function outerDifference_r
        
        pure function outerDifference_i(a, b)
                integer, intent(in) :: a(:), b(:)
                integer :: outerDifference_i( size(b), size(a) )
!REMOVE                outerDifference_i = spread( a, 2, size(b) ) &
!                                  - spread( b, 1, size(a) )
                outerDifference_i = spread( a, 1, size(b) ) &
                                  - spread( b, 2, size(a) )
        end function outerDifference_i

        pure function outerProduct_4(a, b)
                real(4), intent(in) :: a(:), b(:)
                real(4) :: outerProduct_4( size(b), size(a) )
                outerProduct_4 = spread( a, 1, size(b) ) &
                                  * spread( b, 2, size(a) )
        end function outerProduct_4
        
        pure function outerProduct_8(a, b)
                real(8), intent(in) :: a(:), b(:)
                real(8) :: outerProduct_8( size(b), size(a) )
                outerProduct_8 = spread( a, 1, size(b) ) &
                                  * spread( b, 2, size(a) )
        end function outerProduct_8
        
        pure function outerProduct_i(a, b)
                integer, intent(in) :: a(:), b(:)
                integer :: outerProduct_i( size(b), size(a) )
!REMOVE                outerProduct_i = spread( a, 2, size(b) ) &
!                                  * spread( b, 1, size(a) )
                outerProduct_i = spread( a, 1, size(b) ) &
                                  * spread( b, 2, size(a) )
        end function outerProduct_i

        ! Returns upper triangular array mask:
        pure function upperTriangular_0(m, n)
                integer, intent(in) :: m, n
                logical :: upperTriangular_0(m, n)

                upperTriangular_0 = (outerDifference( &
                        arithmeticVector(1, 1, n), &
!REMOVE                        arithmeticVector(1, 1, n) ) >= 0)
                        arithmeticVector(1, 1, m) ) <= 0)
        end function upperTriangular_0

        pure function upperTriangular_up(m, n, up)
                integer, intent(in) :: m, n
                logical :: upperTriangular_up(m, n)
                ! moves diagonal # of steps up:
                integer, intent(in) :: up

                upperTriangular_up = (outerDifference( &
                        arithmeticVector(1, 1, n), &
!REMOVE                        arithmeticVector(1, 1, n) ) >= up)
                        arithmeticVector(1, 1, m) ) <= -up)
        end function upperTriangular_up

        ! Returns lower triangular array mask:
        pure function lowerTriangular_0(m, n)
                integer, intent(in) :: m, n
                logical :: lowerTriangular_0(m, n)

                lowerTriangular_0 = (outerDifference( &
                        arithmeticVector(1, 1, n), &
!REMOVE                        arithmeticVector(1, 1, n) ) <= 0)
                        arithmeticVector(1, 1, m) ) >= 0)
        end function lowerTriangular_0

        pure function lowerTriangular_up(m, n, up)
                integer, intent(in) :: m, n
                logical :: lowerTriangular_up(m, n)
                ! moves diagonal # of steps up:
                integer, intent(in) :: up

                lowerTriangular_up = (outerDifference( &
                        arithmeticVector(1, 1, n), &
!                        arithmeticVector(1, 1, n) ) <= up)
                        arithmeticVector(1, 1, m) ) >= -up)
        end function lowerTriangular_up

        pure subroutine copyToLowerTriangleDeprecated(squareMatrix)
                real(D), intent(inout) :: squareMatrix(:, :)

                integer :: i

                do i = 1, size(squareMatrix, 1) - 1
                        squareMatrix( i, i + 1:size(squareMatrix, 2) ) = squareMatrix( i + 1:size(squareMatrix, 1), i)
                end do
        end subroutine copyToLowerTriangleDeprecated

        pure subroutine copyToUpperTriangleDeprecated(squareMatrix)
                real(D), intent(inout) :: squareMatrix(:, :)

                integer :: i

                do i = 1, size(squareMatrix, 1) - 1
                        squareMatrix( i + 1:size(squareMatrix, 2), i ) = squareMatrix( i, i + 1:size(squareMatrix, 1))
                end do
        end subroutine copyToUpperTriangleDeprecated

        integer pure function iminloc_4(a)
                real(4), intent(in) :: a(:)

                iminloc_4 = sum( minloc(a) )
        end function iminloc_4

        integer pure function iminloc_8_1D(a)
                real(8), intent(in) :: a(:)

                iminloc_8_1D = sum( minloc(a) )

        end function iminloc_8_1D
 
        pure function iminloc_8_2D(a)
                real(8), intent(in) :: a(:, :)
                integer :: iminloc_8_2D( size(a, 2) )
                iminloc_8_2D = sum( minloc(a) )
        end function iminloc_8_2D
 
!        integer pure function iminloc_16(a)
!                real(16), intent(in) :: a(:)
!                iminloc_16 = sum( minloc(a) )
!        end function iminloc_16

        integer pure function imaxloc_4(a)
                real(4), intent(in) :: a(:)
                imaxloc_4 = sum( maxloc(a))
        end function imaxloc_4

        integer pure function imaxloc_8(a)
                real(8), intent(in) :: a(:)
                imaxloc_8 = sum( maxloc(a) )
        end function imaxloc_8

!        integer pure function imaxloc_16(a)
!                real(16), intent(in) :: a(:)
!                imaxloc_16 = sum( maxloc(a) )
!        end function imaxloc_16

        logical function even(n)
                integer, intent(in) :: n
!                even = REAL(n / 2) == ( REAL(n) / 2.0 )
                even = (n / 2) * 2 == n
        end function even

        ! Returns true if all elements of v are true:
        logical function true(v)
                logical, intent(in) :: v(:)
                integer :: i
                true = .true.
                do i = 1, size(v)
                        true = true .and. v(i)
                end do
        end function true

        logical function odd(n)
                integer, intent(in) :: n
!                odd = REAL(n / 2) /= ( REAL(n) / 2.0 )
                odd = (n / 2) * 2 /= n
        end function odd

        subroutine allocate_r_1D_ptr(ptr, size, name)
                real(D), pointer :: ptr(:)
                integer, intent(in) :: size
                character(len = *), intent(in) :: name
                integer :: ierr

                if( associated(ptr) ) then
                        write(0, *) 'ERROR:', name, 'is already allocated.'
                        stop
                end if

                allocate( ptr(size), stat = ierr )
                if(ierr /= 0) then
                        write(0, *) 'ERROR: Allocation of', name, 'failed.'
                        stop
                end if
        end subroutine allocate_r_1D_ptr

        subroutine checkAlloc(ierr, name)
                integer, intent(in) :: ierr
                character(len = *) :: name
                if(ierr /= 0) then
                        write(0, *) 'ERROR: Allocation of', name, 'failed.'
                        stop
                end if
        end subroutine checkAlloc

        subroutine deallocate_r_1D_ptr(ptr, name)
                real(D), pointer :: ptr(:)
                character(len = *), intent(in) :: name
                integer :: ierr

                if( .not.associated(ptr) ) then
                        write(0, *) 'WARNING: tried to deallocate unallocated', name
                else
                        deallocate( ptr, stat = ierr )
                        if(ierr /= 0) then
                                write(0, *) 'ERROR: Deallocation of', name, 'failed.'
                                stop
                        end if
                end if

        end subroutine deallocate_r_1D_ptr

 ! pgf90 won't take allocatable arrays as dummy arguments
#ifdef IFORT
        subroutine deallocate_array(A, name)
                real(D), allocatable :: A(:)
                character(len = *) :: name
                integer :: ierr
                if( allocated(A) ) then
                        deallocate(A, stat = ierr)
                        if(ierr /= 0) then
                                write(0, *) 'ERROR: Deallocation of', name, &
                                            'failed.'
                                stop
                        end if
                end if
        end subroutine deallocate_array
#endif ! IFORT

        logical pure function almostEqual(val1, val2, threshold)

                real(D), intent(in) :: val1, val2, threshold

                if(val1 /= 0.0_D .and. val2 /= 0.0_D) then
                        almostEqual = (( abs(val1 - val2) / val1 ) <= threshold)
                else
                        almostEqual = ( abs(val1 - val2) <= threshold )
                end if

        end function almostEqual

end module mMisc
